Density profile of interacting Fermions in a one-dimensional optical trap 
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The density distribution of the one-dimensional Hubbard model in a harmonic trapping potential 
is investigated in order to study the effect of the confining trap. Strong superimposed oscillations 
are always present on top of a uniform density cloud, which show universal scaling behavior as a 
function of increasing interactions. An analytical formula is proposed on the basis of bosonization, 
which describes the density oscillations for all interaction strengths. The wavelength of the dominant 
oscillation changes with interaction, which indicates the crossover to a spin-incoherent regime. Using 
the Bethe ansatz the shape of the uniform fermion cloud is analyzed in detail, which can be described 
by a universal scaling form. 
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Ultra-cold gases in optical traps and lattices have be- 
come a promising tool for simulating strongly correlated 
systems with a full control of all relevant parameters 
[ij. While the first simulations were mostly made on 
bosonic setups, ultra-cold fermions are by now also well 
established [2[- In order to simulate interacting electron 
systems such as Hubbard- type models, fermionic atoms 
with two different hyperfine states are used in order to 
represent the two spin channels 0]. It is therefore possi- 
ble to test theoretical predictions even for systems that 
are less common or hard to produce in nature, such as 
perfectly clean isolated one-dimensional (ID) quantum 
wires. However, the experimental setup will always pos- 
sess a smoothly varying potential due to the intensity 
profile of the laser beams, usually forming a harmonic 
confinement. 

Recent experimental developments have made it pos- 
sible to locally probe the density profile of ultra-cold 
atomic condensates directly in space using optical mi- 
croscopy 0] or electron beam scanning In this work 
we therefore want to provide a detailed theoretical quan- 
titative analysis of the ID fermion density profile as a 
function of the interaction strength and the confining po- 
tential, which in turn can be used to analyze interaction 
effects from the experimental signals. 

The fermion density can generally be characterized in 
terms of two distinct features, namely the overall size of 
the cloud on the one hand and superimposed density os- 
cillations on the other hand. From works on quantum 
wires and quantum dots it is well known that density 
oscillations may appear from reflections at sharp edges 
and boundaries, which are due to interference (Friedel 
oscillations) and/or localization (Wigner crystallization) 
[sl-fioj. However, it is not a priori clear how these os- 
cillations are modified if a harmonic potential is present 
as a confinement. In this work we now show that the 
oscillations remain strong in a harmonic trap with inter- 
actions, despite the lack of any sharp edges which may 
cause Friedel oscillations. In a pioneering work [T]| from 
1993, Schulz predicted Akp density correlations to dom- 



inate, which he called a one-dimensional "Wigner crys- 
tal", that only occurs in case that the interaction pa- 
rameter takes on rather extreme values, which cannot be 
reached in a short-ranged Hubbard model even for infi- 
nite U . In contrast to this expectation, we now find a 
surprising crossover towards rather strong Akp "Wigner 
oscillations" in a trap even for intermediate short-range 
interactions. Moreover, both Friedel and Wigner oscil- 
lations can actually be very well analyzed with the help 
of an analytic formula on the basis of a bosonization ap- 
proach. The overall size and shape of the density cloud 
is also analyzed in quantitative detail, which follows a 
universal scaling form. 

We consider the standard ID Hubbard Hamiltonian 
with an external trapping potential 
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in the limit of large particle separations (small densi- 
ties) relative to the lattice spacing. In this limit [l^] the 
Hamiltonian can also be approximated by the continuous 
problem of fermions with contact interactions 
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where we assume a non-magnetic state with fixed parti- 
cle number N = 2N^ ~ ^N^. The lattice spacing and the 
hopping J are the natural units for this problem which 
arc set to unity in what follows. The condition for large 
particle separation corresponds to Noj ^1. It should 
be noted that the opposite limit of small particle sepa- 
rations (i.e. of the order of the lattice spacing) has been 
studied elsewhere and is governed by a transition to a 
Mott-insulator [Tsl - figj . A good qualitative understand- 
ing of the continuous problem has been achieved using 
density functional methods [l^, [2l| . 
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FIG. 1: (Color online) DMRG data for the fermion density 
n{x) in a harmonic trap with A'' = 30 and lo'^ — 4 ■ 10"'' 
(points) compared to the analytical approximation in Eq. (|12p 
(solid lines). 



1.2h 

l.ll 



I V 



U/VNu 25 



50 



FIG. 2: (Color online) The width Lf{U) of the fermion cloud 
in the trap as determined from the local Bethe ansatz density 
as a function of the scaling variable U/^/WiZ. The differ- 
ent symbols (colors) correspond to many different choices of 
U £ [0,20], uj'^ G [10-^16 X 10"^] and' AT G [10,70]. Inset: 
Effective exponent 5 determined from a fit of rip (a;) to the 
local Bethe ansatz density. 



In order to simulate the Hamiltonian in Eq. ([1]) we 
use the numerical density matrix renormalization group 
(DMRG) [12. While the DMRG is best suited for ho- 
mogeneous systems with open boundary conditions, it 
is also possible to implement the algorithm to describe 
inhomogencous traps as long as the actual system size 
in the simulation is much larger than the spread of the 
confined fermions. 

In Fig. [1] typical density distributions from DMRG are 
plotted for different U in a trap with iV = 30 particles 
and = 4 • 10^^, showing a localized fermion cloud 
with superimposed characteristic oscillations of different 
wavelengths. An analytical approximation to the data 
is also shown (solid lines), which will be derived in the 
following using bosonization and Bethe ansatz methods. 

In order to understand the behavior of the den- 
sity let us first consider the non-interacting case. 
In the continuous limit the wave-functions of the 
single-particle oscillator levels are given by hn{x) = 

^f^{^y^^e-'^^"/''H„{y^x), where ff„(a;) denotes 
the n-th Hcrmite polynomial. The ground state density 
distribution at U — can be calculated as the sum over 
the filled Fermi sea of oscillator levels 

N/2-1 

n{x) = 2 \hn{x)\^ (3) 

for a system containing N electrons. Using an expansion 
around the center of the trap, this function can be well 
described by a simple closed formula 23 25] 

(-1)^/2 cos{2kF{x)x) 
n{x) « no a; p z 2772-' W 

for \x\ < Lp, where the density cloud is given by 



with a Thomas- Fermi size ot Lp = \/N/uj for U = 0. 
The result in Eq. ^ resembles the corresponding ex- 
pression for Fricdel oscillations in a ID box which 
also decay proportional to the reciprocal distance from 
the turning points ±Lp. The slowly varying part of the 
density 710(2;), Eq. ([5]), replaces the normally constant 
filling uq. The period of the oscillations is related to 
the filling, so that the wavevector also becomes position 
dependent and is given by the non-local expression 
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which follows from an expansion of the summed up oscil- 
lator wave- functions in Eq. ([3]). Note that the integration 
in Eq. ([7]) is similar to the usual WKB approximation, 
where the local momentum k{x) w 7rri(x)/2 is integrated 
in space in order to predict the behavior in a changing 
potential. For constant filling the expression in Eq. ([7]) 
reduces to the usual relation fcp = §no. 

Before we analyze the oscillations in the presence of 
interactions, let us first consider the overall shape of 
the density cloud no(a;) for repulsive interactions U > 
0. For a translationally invariant ID Hubbard model, 
the density is known as a function of U and chemical 
potential fi from the non-trivial solution of the Bethe 
ansatz equations (26| . This is in sharp contrast to 
bosonic systems where the density is well described by 
a mean field approach, which can even be applied lo- 
cally to non-uniform systems with the help of the Gross- 
Pitaevskii equation. One promising approach for the 
non-homogeneous fermion system may be to use the ex- 
act solution. In particular, if the external potential is 
slowly varying, the Bethe ansatz density for the chemical 
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FIG. 3: (Color online) Density oscillations around the slowly 
varying part 710(2;) for A'' = 30 and oj^ = 4 - 10~^ compared to 
the analytical form in Eq. (|12|) . 



potential [lix] = /io + w^a;^ could be a good approxima- 
tion for each location x in the trap. Since strong long- 
range correlations exist, it is not a 'priori clear if such 
a local Bethe ansatz approximation with a translation- 
ally invariant system at each point is appropriate, but it 
agrees very well with our DMRG data for all U . This 
is especially surprising near the edges where the filling is 
low and also the discrete energy spectrum should play a 
role. At [/ = this approach corresponds to the den- 
sity in Eq. At [/ — > 00 the interaction induces a 
Pauli principle between spin-up and down electrons, so 
that the system contains effectively twice N = -|- 
non-interacting spin-incoherent particles, which results 
again in the density in Eq. ([5]), but with Lp multiplied 
by a factor of \f2. For intermediate U the total size of 
the cloud therefore becomes interaction dependent with 
^jNjijj < Lf{U) < y/2N/u}. We find that the interac- 
tion dependence of the effective size L p (U) in fact follows 
a universal scaling behavior as a function of U/ \/Nlo as 
shown in Fig. [21 which is related to the scaling behavior 
of the Bethe ansatz equations with U/hq at low filling 
no 0, • The Bethe ansatz density at intermediate U 
does not follow exactly the simple expression in Eq. ([5]). 
but can be approximated by the normalized shape if it is 
raised by a small exponent of (5 < 1.3 as shown in the inset 
in Fig. [2j It should be noted that Fig. [2] gives a quanti- 
tative estimate of the screening cloud for all interaction 
strengths U and trap parameters as long as Noj <C 1. 

We can now subtract the Bethe ansatz estimate for 
the slowly- varying part of the density no(a;) from the 
DMRG data in order to analyze the oscillations as shown 
in Fig. [3l For weak interactions the Friedel-type oscilla- 
tions in Eq. (U) can clearly be seen. At intermediate 
interactions [/ = 4 two dominant wavevectors can be 
observed. At larger U the faster oscillations dominate, 
corresponding to exactly one density maximum per par- 
ticle, which is one typical signature of Wigner crystal 



oscillations. The oscillating signal is quite sensitive to 
the estimate of the uniform density, which has been sub- 
tracted. However, as can be seen in Fig.[3]the oscillations 
are symmetric in the entire trap without any visible bias 
towards positive or negative values, which shows that the 
local Bethe ansatz estimate works well. 

A natural tool for calculating the density oscillations 
with the help of correlation functions in one-dimensional 
systems is bosonization Q . In the presence of a trapping 
potential, bosonization has been considered for interact- 
ing spinless fcrmions before [27l - [30| . For the spinful case 
we will now derive the central definitions of the bosonic 
creation and annihilation operators. Instead of the usual 
left- and right-moving bosons, only one bosonic field each 
for spin and charge (z/ = c, s) is defined 
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where the bosonic annihilation and creation operators 
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are expressed in terms of fermion operators cj. j„ of the 
TO-th oscillator mode that are extended to include non- 
physical states m < (with ± corresponding to v = 
c,s respectively). The number operators iV^ ± are 
canonical conjugate to the zero modes 0^. The auxiliary 
variable u S [— ti", 7r[ should not be confused with the 
position X. Following the usual steps of the bosonization 
procedure [3l|, it is then easy to show that free-particle 
excitations relative to the Fermi edge are reproduced by 
the bosonic Hamiltonian 
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The Fourier transformation of the oscillator levels de- 
fines a canonical auxiliary fermion field which can be 
bosonizcd as a vertex operator as usual 
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The auxiliary field can in turn be used to give a non-local 
expression of the physical fermion fields in terms of the 
bosons 

00 

i^lix) = ^/i„(a::)4,„ 

ri=0 
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This expression can be made approximately local 
in X by noticing that the wave-functions near the 
Fermi edge oscillate roughly as a function of n with 
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FIG. 4: (Color online) Amplitudes Ai,2 (top) and exponents 
cn,2 (bottom) as determined by fitting Eq. (fT2)) to DMRG 
data for different choices of U, to, and A'^. 
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Therefore, ipl-{x) ^ 

arccos(x/LF), which also leads to densities in terms of 
derivatives of the boson field j27l - [30| . Without interac- 
tions this bosonization approximation is in fact more ac- 
curate than for translational invariant systems due to 
the linear oscillator spectrum. However, spinful inter- 
actions become quite complicated in the bosonized lan- 
guage, since general scattering terms appear that are not 
even momentum conserving. 

Although we are not able to solve the system by a 
simple Bogoliubov transformation, the bosonization pic- 
ture is useful since it is reasonable to expect that the 
leading instabilities are again a 2fci? Friedel and a Akp 
Wigner oscillation as for the translational invariant sys- 
tem 0, [HI albeit with a changing wavevector along the 
trap according to Eq. ([7]). We can therefore generalize 
Eq. (HI analogously to the Hubbard model with fixed 



boundary conditions and propose a general ansatz for 
the density in the trap for \x\ < Lp 



n{x) = uq^x) — A\ 
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where Lp{U) is given in Fig. [5] and k-p{x) is given in 
Eq. (O as a function of Lp. The amplitudes Ai^2 and the 
exponents ai_2 are unknown and have to be determined 
from fitting the DMRG data. As can be seen in Fig. [3] 
this formula fits the data extremely well even to within 
the last oscillation near the edge. In Fig. |4] the results 
for the amplitudes and the exponents are shown which 
again follow a scaling law as a function of U /\/ Nui. For 
smaller amplitudes Ai^2 ^ 0.2 the corresponding expo- 
nents in Fig. 2] could no longer be accurately determined. 
A clear increase of the faster Wigner-crystal oscillations 
A2 can be seen with increasing U. At the same time the 
slower Friedel-type oscillations Ai arc suppressed. Both 
decay exponents generally decrease with increasing inter- 
actions which is qualitatively similar to the translational 
invariant Hubbard model |26| . 

In conclusion we have analyzed the detailed behavior 
of the fermion density of the one-dimensional Hubbard 
model in a harmonic trap with the help of bosonization 
and the Bethe ansatz. The proposed analytical formula 
in Eq. ((T2|) and the scaling behavior of the parameters 
in Fig. [21 and Fig. S] provide very accurate predictions for 
the position dependent density n(x) as a function of ar- 
bitrary interaction strengths U and trap parameters in 
the limit Nlu <^ 1. Significant deviations can only be 
observed in the last oscillations near the edge of the den- 
sity cloud. The overall density no(a;) follows a local Bethe 
ansatz approximation and the oscillations in the trap re- 
main strong despite the lack of any hard-wall boundary 
conditions. A crossover from slower Friedel oscillations 
to faster Wigner crystal oscillations can be observed with 
increasing U. We hope that our results will be useful in 
the analysis of future experiments on ultracold fermions 
in a one-dimensional trap with local resolution. At the 
same time the good fit to the proposed analytical for- 
mula in Eq. (|12p strongly suggests that the problem can 
be solved by further analyzing the bosonization formu- 
las in Eqs. ([71 fTT|) in the presence of interactions, which 
might inspire future research on the topic. 
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